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Large-eddy simulation of turbulent flow with 
a surface-mounted two-dimensional obstacle 

By Kyung-Soo Yang 1 AND Joel H. Ferziger 1 


1. Motivation and objectives 

Large-eddy simulation (LES) is an accurate method of simulating complex tur- 
bulent flows in which the large flow structures are computed while small scales are 
modeled. The rationale behind this method is based on two observations: most 
of the turbulent energy is in the large structures, and the small scales are more 
isotropic and universal. Therefore, LES may be more general and less geometry- 
dependent than Reynolds- averaged modeling, although it comes at higher cost. 

Even though LES has been used by many investigators, most research has been 
limited to flows with simple geometry. Here we shall consider a rectangular paral- 
lelopiped mounted on a flat surface. Related flows are those over surfaces protruding 
from submarines (conning towers or control fins), wind flows around buildings, and 
air flows over computer chips, among others. The most distinctive features associ- 
ated with these flows are three dimensionality, flow separation due to protruding 
surfaces, and large scale unsteadiness. As a model flow, we consider a plane channel 
flow in which a two-dimensional obstacle is mounted on one surface (see Fig. 1). 
This relatively simple geometry contains flow separation and reattachment. Flow 
in this geometry has been studied by Tropea & Gackstatter (1985) for low Re and 
Werner & Wengle (1989) and Dimaczek, Kessler, Martinuzzi & Tropea (1989) for 
high Re , among others. 

Recently, Germano, Piomelli, Moin & Cabot (1991) suggested a dynamic subgrid- 
scale model in which the model coefficient is dynamically computed as computation 
progresses rather than input a 'priori. This approach is based on an algebraic iden- 
tity between the subgrid-scale stresses at two different filter levels and the resolved 
turbulent stresses. They applied the model to transitional and fully turbulent chan- 
nel flows and showed that the model contributes nothing in laminar flow and exhibits 
the correct asymptotic behavior in the near-wall region of turbulent flows without 
an ad hoc damping function. This is a significant improvement over conventional 
subgrid-scale modeling. 

Until very recently, use of the dynamic model in complex geometries has been dif- 
ficult owing to the lack of homogeneous directions over which to average the model 
coefficient (see Ghosal et al this volume for a dynamic model applicable to inho- 
mogeneous flows). The present work was accomplished prior to the developments 
of Ghosal et al and accordingly makes use of a combination of time and spatial 
averages in order to determine the model coefficient. The averaging scheme will be 
discussed in more detail in §3. 
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Figure 1. 



In this paper, we perform an LES of turbulent flow in a channel containing a two- 
dimensional obstacle on one wall using a dynamic subgrid-scale model (DSGSM) at 
i?e=3210, based on bulk velocity above the obstacle ( U m ) and obstacle height (/i); 
the wall layers are fully resolved. The low Re enables us to perform a DNS (Case I) 
against which to validate the LES results. The LES with the DSGSM is designated 
Case II. In addition, an LES with the conventional fixed model constant (Case III) 
is conducted to allow identification of improvements due to the DSGSM. 

We also include LES at jRe=82,000 (Case IV) using conventional Smagorinsky 
subgrid-scale model and a wall-layer model. The results will be compared with the 
experiment of Dimaczek et al. (1989). 


2. Formulation 

All variables are nondimensionalized by U m and h. The code uses a nonuniform 
Cartesian staggered grid in a finite-volume approach. The incompressible momen- 
tum equations filtered by a simple volume- average box filter are 


dUj d dp 1 d 2 ^. 

dt dxj' U * dx, Redxjdxj' 


( 1 ) 


where m, u 2 , u 3 (or u, v, w) are velocities in x, (streamwise), x 2 (normal), x 3 
(spanwise) directions (or x, y, z), respectively, and p is pressure. The volume- 
average box filter is defined by 


Ui(x,t ) 
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where ?=(xi , x 2 , x 3 ) and dx'^dx^dx^dx^. The convective term can be rewritten as 


c) d 

— (uTuy) = fcT.iui 11 ] + Lij + Cii + R '})' 


( 3 ) 
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where 


L{j — UiUj UiUj 
Cij = uiu'j + uju'i 

Rij = iijuj. 


( 4 ) 


Li j, Cij , Rij represent Leonard stresses, cross terms, subgrid- scale Reynolds stresses, 
respectively. When a finite-difference scheme of second-order accuracy is used, the 
Leonard stresses are of the same order as the truncation error (Shaanan, Ferziger 
&: Reynolds, 1975). The other terms have to be modeled. 

The governing equations for LES become 

du, , d dP _d_ 1 (Fu, , R s 

dxj ' 3 Re dxjdxj' 


dt 


, dP 
+ “ -dT, 


where 


Tij — Qij ^Qkk&ij 
P = P + \Qkk 
Qij = Rij + Cij. 


( 7 ) 


Here, Sij is the Kronecker symbol. In the present simulation, the eddy viscosity 
model of Smagorinsky (1963) is used: 


= —2vtS 




( 8 ) 


where 



1 > dxii duj v 
2^dxj ^ dxi ' 
= l 2 y0^~ j . 


( 9 ) 


Here, / is a characteristic length scale of the small eddies. In Case III and IV, 
the smaller value of xd and 0.1 A is used for /, where k and d are von Karman’s 
constant and the distance normal to a wall, respectively, and A = ( Ax Ay Az)^ . The 
particular form of in (7) is chosen in order to make both (7) and (8) consistent 
on contraction ( i = j). 

In Case II, l 2 — C s A 2 is dynamically determined following the prescription of 
Germano et al. (1991) as modified by Lilly (1992). When the dynamic model 
is used, C s is an instantaneous and local quantity that can vary wildly in time 
and space. This wide variation results in large negative values of C 3 that lead to 
numerical instability. To avoid this difficulty, averaging is performed in space and 
time. (For an alternate approach see Ghosal et al . , this volume). Spatial averaging is 
done in the homogeneous ( z ) direction first. Then additional averaging is performed 
over nine neighboring grid points with the point for which the averaging is carried 
out at the center, using volume weighting, in order to obtain an averaged value of 
C 9 at a given inner grid point. In the near-wall region, averaging is done only in the 
direction parallel to the solid wall, i.e. using three points. It is necessary to repeat 
this process to smooth C 9 sufficiently. Germano et al (1991) found the optimum 

value of the ratio, A/ A, to be two, a value we adopted. 
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3. Numerical method 

To advance the solution in time, a fractional step method (Kim and Mom, 1985) 
is employed. The time- advancement of the momentum equation is hybrid; the 
convective terms are explicitly advanced by a third-order Runge-Kutta scheme and 
the viscous terms implicitly by Crank-Nicolson method: 


k — 1 

= (a* + + fikHui - uf" 1 ) 


7 *AT(uf- 1 ) - CkN ( «?- a ) - (a k + 0 k ) 
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In the expressions for X and iV’(uj), summation is performed on the index j only. 
The momentum equation is time- advanced without implicit pressure terms and 
then projected onto a divergence-free space by introducing <f> that obeys a Poisson 
equation. The latter is solved by a multigrid method which is very flexible and 
more efficient than a number of competitors. For spatial discretization, second- 
order accurate central differencing was used. All terms in the model except the 
cross derivatives are treated implicitly in all three directions to avoid restrictions 
on time steps. The code is well vectorized; a speed of 150 MFLOPS has been 
achieved on CRAY Y-MP/832. 



LES of an obstacle flow 


101 



FIGURE 2. Schematic drawing of SR regions. 


Case 

Re 

Xr 

X r Y r X F Y f 

i 

3310 

6.42 

1.21 0.35 1.51 0.28 

ii 

3210 

6.80 

1.13 0.36 1.51 0.37 

hi 

3330 

7.01 

1.76 0.28 1.35 0.40 


Table 1. Comparison of various SR lengths. 


4. Results and discussion 

]).l Choice of parameters and boundary conditions 

The values of the geometric parameters in all four cases are hjH— 0.5, W/h= 2, 
and L/h— 1, where H , W, and L are channel height, spanwise width of the obstacle 
and channel, and obstacle streamwise length, respectively (see Fig. 1). The inlet 
and the outlet are located at x=0 and x=31, respectively, and the obstacle is placed 
between x=10 and £—11. 

In Cases II and III the center of the control volume adjacent to the wall was placed 
at Ay=0.0086 from horizontal walls except on the top of the obstacle where the 
nearest center was placed at Ay— 0.0046. The corresponding distances for Case I and 
Case IV axe 0.005 and 0.05 from horizontal walls and 0.0036 and 0.025 from the top 
of the obstacle. On the forward-facing wall, the first grid points are at Ax=0.0045 
for Cases II and III and at 0.0033 and 0.025 for Cases I and IV respectively. On 
the backward-facing wall, the first grid points are at Ax=0.014 for Cases II and 
III and at 0.0045 and 0.025 for Cases I and IV respectively. The grid is densely 
packed around the obstacle and near the channel walls and geometrically stretched 
in the other regions. The number of control volumes in the x, y, and z directions 
are 112 x 48 x 40 for Cases II and III, 272 x 64 x 64 for Case I, and 96 x 32 x 32 for 
Case IV. Grid refinement shows that the spatial resolution is adequate; using more 
control volumes shows improvement, but the difference is insignificant. 

In all cases, periodic boundary conditions were employed in the homogeneous 
(z) direction. At the walls, no-slip boundary conditions were imposed except for 
Case IV where a wall-layer model was employed. We also apply periodic boundary 
conditions in the x direction in order to avoid any uncertainty related to outflow 
boundary condition which has been an area of controversy and to assure a reasonable 
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FIGURE 3(a). Averaged wall shear stress at the lower wall: o, DNS; A , LES with 
DSGSM; +, LES with C,=0.01. 



FIGURE 3(6). Averaged wall shear stress at the upper wall: o, DNS; A , LES with 
DSGSM; +, LES with C,=0.01. 
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FIGURE 4. Averaged velocity profiles at £=12; (a) streamwise (t/), ( b ) normal 
(' V ): o, DNS; a , LES with DSGSM; + , LES with C,=0.01. 

flow approaching the obstacle. Therefore, we are actually simulating an infinitely- 
long channel flow with a periodic array of obstacles. To minimize the interaction 
between “neighboring” obstacles, the long streamwise computational domain (31 h) 
is employed. 

Since the pressure difference between the inlet and the outlet is fixed, Rc is 
slightly different in each case. To match the Reynolds numbers of the various cases 
as closely as possible, we adjusted the pressure difference slightly. The second 
column of Table 1 shows Re for each low-i?e case. The 3% variation in Re should 
be kept in mind in the comparisons below. The high -Re case will also be compared 
with the experiment of Dimaczek et ai (1989) at slightly different Re (Re— 84,000). 
After an initial transient period, the flow becomes fully turbulent and sustained. 
Then, an averaging is performed in the homogeneous direction and in time in order 
to obtain averaged quantities. The time- averaging was taken over 27 characteristic 
time units (h/U m ) for low-Re cases and 38 units for Case IV. 

4 2 Averaged flow field at Re-3210 

The flow develops several separation and reattachment (SR) zones near the ob- 
stacle. Figure 2 shows schematic contours of JJ =0 (U and V represent averaged 





a O 



FIGURE 5. Averaged turbulent fluctuation profiles; (a) streamwise (u' 2 ) at a 1 — 12, 
(6) normal (i/ 2 ) at .r=12, (c) spanwise {w' 2 ) at x=ll: o, DNS; A , LES with 
DSGSM; + , LES with C,=0.01. 
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values of u and v , respectively). The + and - signs indicate regions of positive and 
negative £/, respectively. In addition to the primary SR zones upstream (PSRU) 
and downstream (PSRD) of the obstacle, there are secondary SR zones upstream 
(SSRU, at the front corner), and downstream (SSRD, near the rear corner, bigger 
than SSRU) of the obstacle. A tertiary SR zone is discernible at the downstream 
corner of the obstacle (TSRD). For the given geometry and Re, reattachment does 
not occur on top of the obstacle in any of the three low-Re cases. The reattachment 
length of PSRD is denoted by Xr. The separation length and reattachment length 
of PSRU are represented by Xf and Yp, respectively, and those of SSRD by X r 
and Y r , respectively. Table 1 gives computed values of those lengths in units of h 
for each case. Case I, the DNS, is the most accurate simulation. Its Xr differs by 
roughly 5% from the value of 6.1 determined in the experiment of Tropea & Gack- 
statter (1985). Although the aspect ratio ( Lfh ) of the experimental obstruction is 
somewhat larger than in the simulation, the DNS value of Xr falls within the range 
of experimental error (quoted as 6%). Tropea & Gackstatter did not report other 
SR lengths. Case II is significantly more accurate than Case III for every length 
scale. 

Figure 3 presents nondimensionalized shear stress (r^) on the lower (Fig. 3(a)) 
and upper walls (Fig. 3(b)). In Fig. 3(a), the values of r w between .t= 10 and x=ll 
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FIGURE 7. Regions of instantaneous negative u; solid, positive; dash, negative; 
thick solid, 0; increment, 0.016: (a) t— 0; ( 6 ) (c) t— 2At. 


are for the top surface of the obstacle. The large variations in the values of r w 
on the top of the obstacle reflect the complexity of the flow in that region. The 
t w predicted by Case II agrees better with Case I in PSRD than does Case III, 
especially for 11 < x < 13 and far downstream ( x > 20). Case II also gives better 
results on the upper channel wall (Fig. 3(6)). The large jr^l near x=T0 is caused by 
flow acceleration due to the sudden contraction in flow passage. Better agreement 
for 11 < x < 15 and in the “channel region” (x < 7.5 or x > 25) are also noticeable. 

Profiles of U and V at a selected streamwise location (x— 12, just downstream of 
the obstacle) are shown in Figs. 4(a) and Fig. 4(6), respectively. In both figures, 
the DSGSM gives a significant improvement over the Smagorinsky model in the 
reversed flow region (y < 0.75). 

Profiles of averaged turbulent fluctuations in the streamwise (u /2 ), normal (i> /2 ), 
and spanwise {w f2 ) directions at selected streamwise locations are presented in 
Figs. 5(a), 5(6), and Fig. 5(c), respectively. It should be noted that the LES re- 
sults represent only the fluctuations in the resolved (grid-scale) velocity field. The 
subgrid-scale contribution is small at this low Re, The dynamic model gives an 
overall improvement for u ' and v ! , but not for w . 

Figure 6 shows profiles of C 9 (x, y) at three different streamwise locations (.r=10.8, 
12, 15). Obviously, C s depends upon grid used and the type of averaging in space 
and time. There is a sharp gradient near y=l where the control volumes are densely 
clustered to resolve the flow above the obstacle. Without an arbitrary damping 
function, C 3 vanishes at the walls and even takes some small negative values near 
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FIGURE 8. Regions of instantaneous negative u; solid, positive; dash, negative; 
thick solid, 0; increment, 0.016: (a) t=0; ( b ) t=At\ (c) t=2At. 

the upper wall. 

4.3 Instantaneous flow field at Re =3210 
4*3.1 Reversed flow regions 

Figure 7 shows contours of u are presented at one x-y plane at three different 
times with a time interval of At =1.61. For convenience, the time for Fig. 7(a) is 
designated as t= 0 , and subsequent figures will be referred relative to that time. 
Figures 7(a)— (c) show how unsteady the flow is. Near the mean reattachment 
point of PSRD (6.8/1 downstream of the obstacle), u is small and oscillating in sign. 
Intense unsteady free-shear layers formed downstream of the obstacle are noticeable. 
Intermittent separation on the upper channel wall is observed near the streamwise 
location of mean reattachment (Fig. 7(c)). 

Figures 8(a)-(c) show contours of u at the first grid point away from the lower 
channel wall at three different times. The instantaneous separation and reattach- 
ment lines are far from two-dimensional although the obstacle is geometrically two- 
dimensional. Secondary and tertiary flow regions are present near the obstacle at 
this Re and are highly unsteady. 

Particle trace studies were also performed. A videotape displaying this data is 
available by request to the authors. 

4.4 LES at Re =82,000 

The wall-layer model we used is a variation of one proposed by Ciofalo and Collins 
(1989) for k-e modeling of turbulent recirculating flows. It retains the form of the 
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Figure 9. Streamwise velocity profiles: (a) x=9.6, ( b ) ;r=10.8, (c) x=15: o, exp. 
A , 96x32x32; +, 128x48x40. 
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wall function but allows the nondimensional thickness of the viscous sublayer to be 
a function of the local turbulence intensity. 

Figure 9 presents streamwise velocity profiles normalized by 2 £/ m and averaged 
in 2 and t at three streamwise locations. Velocity profiles from coarse (96 x 32 x 32) 
and fine (128 x 48 x 40) grid simulations are shown along with the experimental 
one. The profiles are relatively well resolved. Only a slight improvement is obtained 
on the fine grid. The discrepancy near the top surface of the obstacle (Fig. 9(b)) 
is believed to be due to the wall-layer model. Averaged and normalized velocity 
fluctuations at selected x locations are shown in Fig. 10. Numerical results predicts 
higher values in the high speed regions near the upper wall. 

5. Summary 

A large-eddy simulation of low-Reynolds-number turbulent flow in a channel with 
a two-dimensional obstacle on one wall was presented with the wall layers fully re- 
solved. The subgrid-scale model coefficient was computed dynamically. The results 
obtained were compared with a DNS and showed that the dynamic model yields 
better results than conventional LES with a fixed model constant. This demon- 
strates the value of the dynamic subgrid- scale model for computing complex flows. 
A high-Reynolds-number LES using a conventional Smagorinsky model with a fixed 
model constant was also included. The results are consistent with the experiment of 
Dimaczek et al. (1989). Application of the dynamic subgrid-scale model to high-Re 
flows is currently under investigation. 
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